function zero = sdp_solve(in,par)

wbar = in(1);
chi  = in(2);

% Unpack parameters
mu_ = par.mu_;
model = par.model;
mfreq = par.mfreq;
estim = par.estim;
show = par.show;

freq_target = 0.184;
size_target = 0.1051;
kurt_target = 4.1;
N_target    = 1/3;

menucost=NaN; lbar=NaN; % initialize parameters
ksi=NaN; noise=NaN; %#ok<*NASGU> % initialize legacy code params

switch model
        case 2  % CALVO 
                lbar = 0.13;
                if estim
                   std_eA = in(2);
                else
                   std_eA = 0.002298088835207;
                end 
        case 3  % GOLOSOV-LUCAS 
                if estim
                    menucost =  in(2); % fixed menu cost 
                    std_eA = in(3);
                else
                    menucost = 0.013818996429905;
                    std_eA = 0.001087678420834;
                end
        case 4  % Woodford 
                lbar = freq_target;        %set to the target frequency
                if estim
                    menucost =  in(2); % fixed menu cost 
                    std_eA = in(3);
                    alpha  = in(4);   %info cost
                else
                    menucost = 0.480356448604763; %chi6: 0.2417;   %eta033:   0.2417;  %0.5;   %eta0: 
                    std_eA = 0.060650925776505;  %chi6: 0.0606;  %eta033:  0.0606;   %0.06;   %eta0: 0.0604     
                    alpha  = 5.383174036965086;  %%chi6: 2.7252;  %eta033:  2.7252;   %5;   %eta0: 0.9579  % info cost (alpha)  
                end                
end
    
% Macro parameters
beta  = (1.03)^(-1/mfreq); % time discount factor
gamma = 2;                 % CRRA coefficient
%chi   = 1.5;  %6            % labor supply parameter
theta = 6;                 % Dixit-Stiglitz elast. of subst.
nu    = 1;                 % log money demand parameter
eta   = 1/3; %0/3;               % degree of real rigidity (decreasing returns)

% Autocorrelations 
rho  = 1;   % Idiosyncratic productivity process  
phiR = 0.9; % Money growth shock
phiA = 0;   % TFP shock

if lbar>1 || lbar<0 || rho<0 || rho>1 || menucost<0 || std_eA < 0
    zero = NaN(size(in)); % tells solver to check another point
    disp('Limits reached on parameters')
    return;
end

nump = 101;   %51 %1001; % Number of price grid points

numstddevs = 9; %9 6;  % Number of std devs around midpoint for idiosyncr. prod.

[T2, Pgrid] = TRP_AR(rho, 0, std_eA^2, nump, numstddevs);

PMAT = exp(Pgrid);  % price grid in levels

pstep = Pgrid(2)-Pgrid(1);  % distance between price grid points

mu = (1 + mu_/100).^(1/mfreq); % gross inflation at model frequency

R = Rmatrix(nump, mu, pstep);  

Cbar = (wbar/chi)^(1/gamma);  
%PAYOFFMAT = Cbar.*PMAT.^(-theta).*(PMAT-wbar);

%PDist1 = zeros(nump,1);
%PDist1((nump+1)/2) = 1;
%PD1 = sum(sum((PMAT.^(-theta)).*PDist1));
%Ybar = Cbar/(1-(eta/(1-eta)*wbar)^(1-eta)*PD1);
%Yi = PMAT.^(-theta)*Ybar;

PAYOFFMAT = Cbar*(PMAT.^(1-theta) - wbar*PMAT.^(-theta/(1-eta)));

%PAYOFFMAT = Yi.*(PMAT - wbar); 

v_iter  

if exitflag==1
    
    Pdist=zeros(nump,1);    
    
    Pdist(ceil(nump/2),:) = 1;
   
    Pdist=Pdist/(sum(sum(Pdist)));   
   
    PdistDIFF=inf;           % reset difference of Pdist
    Piter=0;                 % reset counter
    while (PdistDIFF>sqrt(eps) && Piter<1000)  % iterate to convergence 
        Piter=Piter+1;                    % increment counter
        Pdist_sh = R*T2*Pdist;            % distribution after shocks
        PdistNew = (1-lambda).*Pdist_sh + P.*(ones(nump,nump)*(lambda.*Pdist_sh)); 
        PdistDIFF=nump*max(max(abs(PdistNew-Pdist)));   % sup norm 
        Pdist=PdistNew;                                      % updating Pdist   
    end
    if PdistDIFF>sqrt(eps)
       disp('Pdist not converged!')
       sound(sin(1:300))
    end
    Pdist=Pdist/(sum(sum(Pdist)));   

    calcstats;
    
    if par.comp 
        zero(1) = wbar;
        zero(2) = freqpchanges;
        zero(3) = PD;
        zero(4) = U;   
    elseif show 
        printstats
        plotfigs
        zero = (1-sum(sum(PMAT.^(1-theta).*Pdist))^(1/(1-theta)));
    else
        zero(1) = (1-sum(sum(PMAT.^(1-theta).*Pdist))^(1/(1-theta))); 
        zero(2) = N_target^(eta)-Cbar*PD;
        if estim
           if model == 2
               zero(3) = MeanAbsPchange-0.09;
           elseif model == 3
               zero(3) = freqpchanges-0.13;
               zero(4) = MeanAbsPchange-0.09;
           elseif model == 4
               zero(3) = log(freqpchanges/freq_target);
               zero(4) = log(MeanAbsPchange/size_target);
               zero(5) = log(KurtPchanges/kurt_target);
           end
        end
    end

    if  par.dodyn
        
        dynsolveklein
        
        TT = 60;             
   
        time1Rshock = jacstep; % SET MONEY IMPULSE
        time1TFPshock = 0;
        
        % SPECIFY MONEY SHOCK PROCESS for periods 1:TT
        Rshocks   = [time1Rshock   zeros(1,TT-1)];  % shock + zeros
        TFPshocks = [time1TFPshock   zeros(1,TT-1)];  % shock + zeros
        scalefactor = abs(1/Rshocks(1));    
    
        distsim
        computeirf; 
        plotirf
      
    end

end

end